load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
  fread   = addfile("data/pr_timeseries_1951-2020_JJAS_GPCC.nc","r")
  prN     = fread->prN; original precipitation averaged over northern HMA
  prN     = runave_n_Wrap(prN,9,0,0)
  prS     = fread->prS; original precipitation averaged over southern HMA
  prS     = runave_n_Wrap(prS,9,0,0)

  fread    = addfile("data/pr_eof_1951-2020_JJAS_GPCC.nc","r")
  ts       = fread->ts(0:1,:)
  eof      = fread->eof_regress(0:1,:,:) 

  rad      = 4.*atan(1.0)/180.0 
  clat     = eof&lat
  clat     = cos(clat*rad)
  clat!0   = "lat"
  clat&lat = eof&lat

  pre1  = wgt_areaave_Wrap(eof(0,{32:40},{80:102}),clat({32:40}),1.0,0)
  pre2  = wgt_areaave_Wrap(eof(1,{32:40},{80:102}),clat({32:40}),1.0,0)
  dataN = pre1*ts(0,:)+pre2*ts(1,:); reconstructed precipitation averaged over nothern HMA

  pre1 := wgt_areaave_Wrap(eof(0,{26.5:32},{85:102}),clat({26.5:32}),1.0,0)
  pre2 := wgt_areaave_Wrap(eof(1,{26.5:32},{85:102}),clat({26.5:32}),1.0,0)
  dataS = pre1*ts(0,:)+pre2*ts(1,:); reconstructed precipitation averaged over southern HMA

  print(escorc(prN,dataN))
  print(escorc(prS,dataS)); correlation coefficients


end   